library(INLA)
library(glmmBUGS)
library(R2WinBUGS)

doEverything=T # compute data and do MCMC

NindivPerCommunity = 5
Ncomm = 100
regionsd = 0.2
muNatural = 1/9
mu = log(muNatural/(1-muNatural))


if (doEverything) {

result = data.frame(
  group = 1:Ncomm,
  regionRandomEff = rnorm(Ncomm, 0, regionsd),
  N=NindivPerCommunity
  )
  


result$probsLogit = mu + result$regionRandomEff 
             
result$probs = exp(result$probsLogit) / (1+exp(result$probsLogit))             

result$smoke = rbinom(Ncomm, NindivPerCommunity, result$probs)

subdata = result

forBugs = glmmBUGS(smoke + N ~1, effects="group", family="binomial", data=subdata)

forBugs$startingValues$vars$group = 0.3
forBugs$startingValues$Rgroup = rnorm(length(forBugs$startingValues$Rgroup),
   forBugs$startingValues$intercept, 0.1)

startingValues = forBugs$startingValues

source("getInits.R")

fromBugs = bugs(forBugs$ragged, getInits,
   parameters.to.save = names(getInits()),
   model="model.bug",
   n.chain=3, n.iter=100000,
   n.burnin=100, n.thin=10, debug=T,
   working.directory=getwd())
  
fromBugs$sims.matrix = fromBugs$sims.list = NULL
save(subdata, forBugs, fromBugs, file="bugsResult.RData")
} else {
load ("bugsResult.RData")
}

bugsResult = summaryChain(restoreParams(fromBugs, forBugs$ragged))

bugsResult = rbind(bugsResult$scalars, bugsResult$betas)

if (F){

forBugs = list(y = result$smoke, N=result$N, Ngroup=Ncomm)

getInits = function() {
list(intercept=0, SDgroup = 1, Rgroup =rep(0, Ncomm))
}

fromBugs = bugs(forBugs, getInits, 
  parameters.to.save = c("Tgroup",names(getInits())),
model="modelInla.bug",
  n.chain=3, n.iter=50000, 
  n.burnin=100, n.thin=20,   
  working.directory=getwd() )
fromBugs$sims.matrix = fromBugs$sims.list = NULL  
 save(result, forBugs, fromBugs, file="bugsResult.RData")
}

<<<<<<< .mine
pdf(inlabug)
=======
pdf("inlabug")
>>>>>>> .r149
Sprec =c(  .1, 0.001, .001)
Sprec2 = c(.1, .1, 0.001)
resInla =list()
for(D in 1:length(Sprec)) {
  formula=smoke~f(group, model="iid", param=c( Sprec[D], Sprec2[D]))

  resInla[[D]] = 
	inla(formula,data=result,family="binomial",Ntrials=N,
	quantiles = c(0.025,0.975))

}

par(mfrow=c(2,1))

 hist(fromBugs$sims.array[,,"intercept"], breaks=100, prob=T)
for(D in 1:length(Sprec)) {
 lines(resInla[[D]]$marginals.fixed$intercept, col=D)              
}
 abline(v=mu)           
legend("topleft", legend=as.character(Sprec), 
	lty=rep(1,length(Sprec)), 
	col=seq(1, length(Sprec)))         

     


 hist(fromBugs$sims.array[,,"Tgroup"], breaks=100, prob=T)
for(D in 1:length(Sprec)) {
 lines(resInla[[D]]$marginals.hyperpar[[1]], col=D)              
}
 abline(v=regionsd^(-2))           
 dev.off()
 
 